TK!2
=v
#1
:n

:s
b
:c
*** The Quartic Equation Solver ***
#2
:n
a
:v
1
:s
i
#3
:n
b
:v
-10
:s
i
:c
Coefficients of quartic equation
#4
:n
c
:v
35
:s
i
:c
in general form:
#5
:n
d
:v
-50
:s
i
:c
a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
#6
:n
e
:v
24
:s
i
#7
:n
x1
:v
3.999999999999997
:s
o
#8
:n
x2
:v
3.000000000000012
:s
o
:c
real roots of
#9
:n
x3
:v
1.999999999999988
:s
o
:c
quartic equation
#10
:n
x4
:v
1.000000000000003
:s
o
:c
 
#11
:n
Re1
:s
b
#12
:n
Im1
:s
b
:c
complex roots of quartic equation
#13
:n
Re2
:s
b
:c
  (one of each conjugate pair)
#14
:n
Im2
:s
b
#15
:n
îx1
:v
-1.598721155460225E-14
:s
o
:c
    e      t
#16
:n
îx2
:v
-2.309263891220326E-14
:s
o
:c
    r      e
#17
:n
îx3
:v
-2.309263891220326E-14
:s
o
:c
    r      r
#18
:n
îx4
:v
-1.598721155460225E-14
:s
o
:c
    o      m
#19
:n
îRe1
:s
b
:c
    r      s
#20
:n
îIm1
:s
b
#21
:n
îRe2
:s
b
#22
:n
îIm2
:s
b
#23
:n
y
:s
b
:c
Used for evaluation 
:a
y
#24
:n
x
:v
0
:s
i
:c
of  y = f(x)
:a
x
#27
:n

:s
b
:c
>> See Rule Sheet for Notes <<
=u
=l
#1
:n
x
:v
0
.2
.4
.6000000000000001
.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
3.2
3.4
3.6
3.8
4
4.2
4.4
4.600000000000001
4.800000000000001
5
E
#2
:n
y
:v
24
15.3216
8.9856
4.569599999999999
1.6896
0
-.8064000000000004
-.9984
-.8063999999999999
-.4223999999999999
0
.3456000000000002
.5376000000000002
.5376
.3455999999999996
0
-.4224000000000004
-.8064000000000006
-.9984000000000001
-.8063999999999993
0
1.689600000000002
4.569600000000007
8.985600000000014
15.32160000000003
24
E
=f
#1
:n
quartic
:c
Quartic equation solver
:t
Procedure
:i
a,b,c,d,e
:o
x1,x2,x3,x4,Re1,Im1,Re2,Im2
#2
:n
cubic4
:c
Cubic equation solver for proc. quartic
:t
Procedure
:i
r,s,t
:o
x1
#3
:n
qcheck
:c
Error checker
:t
Rule
:p
a,b,c,d,e
:a
x1,x2,x3,x4,R1,I1,R2,I2
:v
îx1,îx2,îx3,îx4,îR1,îI1,îR2,îI2
#1
:n
quartic
:z
" This function solves 4th order polynomial (quartic) equations

"                4       3       2    
"             a x   + b x   + c x   + d x  +  e   =   0

" for real and complex roots.

" Supporting function:  cubic4

A:= b/a
B:= c/a
C:= d/a
D:= e/a
A2:= A/2
A22:= A2^2
y:= cubic4(-B,A*C-4*D,D*(4*B-A^2)-C^2)   " resolvent cubic equation
R:= sqrt(abs(A22 - B + y))
if abs(y)<1e-7 then y:= 0
if abs(R)<1e-7 then R:= 0
S:= 3*A22 - R^2 - 2*B 
if R=0 then T:= 2*sqrt(abs(y^2-4*D)) else T:= (A*B - 2*C - A*A22)/R
     ST:= S + T
     U:= sqrt(abs(ST))
     if ST<0 goto C12
       x1:= (U - A2 + R)/2
       x2:= x1 - U
       goto 34
C12: Re1:= (R-A2)/2
     Im1:= U/2
34:  U:= sqrt(abs(S-T))
     if S<T goto C34
       x3:= (U - A2 - R)/2
       x4:= x3 - U
       return
C34: Re2:= -(A2+R)/2
     Im2:= U/2
/E
E
#2
:n
cubic4
:z
" This procedure solves a cubic equation in normal form (x^3+r*x^2+s*x+t=0)
" for the first root; it uses Vieta's or Cardano's formulas.
" To solve cubic equation in general form for all three roots, real or
" complex, see TK Library model CUBIC.

p:=s-r^2/3
q:=2*r^3/27-r*s/3 + t
R:=(p/3)^3+(q/2)^2
if R>1e-15 goto Cardano

Vieta:   w:=sqrt(abs(p/3))
         cosine:= round(min((q/2)/w^3, 1), 1e-15)
         phi:=acos(cosine)
         h:= w*cos(phi/3) - r/3
         g:= sqrt(3)*w*sin(phi/3)
         x1:= h + g
         return

Cardano: u3:=-q/2 + sqrt(R)
         v3:=-u3-q
         u:=sgn(u3)*abs(u3)^(1/3)
         v:=sgn(v3)*abs(v3)^(1/3)
         x1:=u+v-r/3
/E
E
#3
:n
qcheck
:z
" This functions may be used for checking the accuracy of roots of 4th
" order polynomial equations evaluated by function  quartic .  The error
" terms  îx1, îx2, ...  differ from zero because of round-off errors.

îx1 = poly(x1,a,b,c,d,e)
îx2 = poly(x2,a,b,c,d,e)
îx3 = poly(x3,a,b,c,d,e)
îx4 = poly(x4,a,b,c,d,e)
(îR1,îI1)=((((a*R1,I1)+(b,0))*(R1,I1)+(c,0))*(R1,I1)+(d,0))*(R1,I1)+(e,0)
(îR2,îI2)=((((a*R2,I2)+(b,0))*(R2,I2)+(c,0))*(R2,I2)+(d,0))*(R2,I2)+(e,0)
/E
E
=r
#1
:r
call quartic(a,b,c,d,e;x1,x2,x3,x4,Re1,Im1,Re2,Im2)
#2
:r
call qcheck(x1,x2,x3,x4,Re1,Im1,Re2,Im2;îx1,îx2,îx3,îx4,îRe1,îIm1,îRe2,îIm2)
#4
:r
y = poly(x,a,b,c,d,e)    " for evaluation of  y = f(x)
:s
C
#6
:r
" Notes:
:s
C
#8
:r
" 1. The real and/or complex roots of a quartic equation are evaluated
:s
C
#9
:r
"    by the procedure 'quartic' which, in turn, calls the a special cubic
:s
C
#10
:r
"    equation solver 'cubic4'.  The direct approach presented here is 
:s
C
#11
:r
"    much simpler and faster then using the Iterative Solver.
:s
C
#13
:r
" 2. There is a slight danger of floating point error propagation through
:s
C
#14
:r
"    the procedure 'quartic', and consequent corruption of the results.
:s
C
#15
:r
"    The rule function 'qcheck' may be used for checking the error terms
:s
C
#16
:r
"    evaluated by applying Horner's rule to given coefficients and calcu-
:s
C
#17
:r
"    lated root values.
:s
C
#19
:r
" 3. The list solving was used and the plot was set up for illustration
:s
C
#20
:r
"    only (press F7).  (Different data may require changes in the  x  list).
:s
C
#22
:r
" 4. See the models CUBIC or BAIRSTOW with TK Library functions for solving,
:s
C
#23
:r
"    respectively, general cubic or higher order polynomial equations
:s
C
#24
:r
"    with real coefficients.
:s
C
=p
#1
:t
Line chart
:d
EGA
:i
""
:s
Yes
:z
X-axis
:g
No
:l
None
:h
""
:v
""
:b
-2.5,10
:x
x
:u
#1
:y
y
:f
1
:p
Lines
:c
"*"
E
=t
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tv,f,0,5,0,1,0
